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Abstract 

A theoretical description for the radial density profile of a finite number of identical charged 
particles confined in a harmonic trap is developed for application over a wide range of Coulomb 
coupling (or, equivalently, temperatures) and particle numbers. A simple mean field approxima- 
tion neglecting correlations yields a density profile which is monotonically decreasing with radius 
for all temperatures, in contrast to molecular dynamics simulations and experiments showing shell 
structure at lower temperatures. A more complete theoretical description including charge correla- 
tions is developed here by an extension of the hypernetted chain approximation, developed for bulk 
fluids, to the confined charges. The results reproduce all of the qualitative features observed in 
molecular dynamics simulations and experiments. These predictions are then tested quantitatively 
by comparison with new benchmark Monte Carlo simulations. Quantitative accuracy of the theory 
is obtained for the selected conditions by correcting the hypernetted chain approximation with a 
representation for the associated bridge functions. 

PACS numbers: 52.27.Lw,52.27.Gr,52.65.Yy 
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I. INTRODUCTION 



Strong correlation effects in ensembles of charged particles are important in many fields 
of physics, including plasmas, the electron gas in solids, and electron-hole plasmas, e.g. 

m 

and references therein. In recent years charged particles spatially confined in trapping poten- 
;ials have attracted considerable interest. Examples are ultracold ions 0, 4|, dusty plasmas 
R y, and electrons in quantum dots |7|. Related studies of expanding neutral plasmas 

riQ 

have been reported as well [8|, 19|. Recently, the ground state structure of Coulomb crystals 
formed by classical charges in a harmonic trap has been studied in detail via molecular dy- 
namics (MD) simulation for Coulomb interaction in [ll| and for Yukawa interaction in 
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12j. The same system is studied here for its fluid phase at finite temperatures. For most 



experimental situations, e.g. in ultracold ions and dusty plasmas, finite temperature effects 
are expected. The objective here is to provide a theoretical description for the structure of 
this system over the full range of temperatures for the fluid phase. 

Among the primary properties of interest is the radial density profile of the particles in 
the trap. A theoretical description of the ground state has been obtained via an energy 
minimization for both Coulomb and Yukawa systems without correlations 13|, and with 
correlations in a local density approximation IJ] . Good agreement with simulations for the 
spatial average density profile was obtained in this way, but a more complete representation 
of correlations is required for the detailed spatial modulation (shells) expected at all except 
the highest temperatures. This is confirmed by the theoretical analysis developed here, 
following This theory is parameterized by the number of particles in the trap, N, 

and the plasma coupling constant, F, measuring the strength of the Coulomb energy for a 
pair of charges in the trap relative to the kinetic energy (defined below). Alternatively, the 
average kinetic energy relative to trap energy or to Coulomb energy can be used to define 
a dimensionless temperature T* = 1/F. The range of values explored is 1 < < 300 and 
< F < 40 (or 0.025 < T* < oo). 

The primary results described here are: 1) an extension of the HNC approximation 
developed for bulk fluids to systems with localized densities, within the context of density 
functional theory; 2) a confirmation that mean field theory (no correlations) predicts a 
structureless, monotonically decreasing radial profile approaching a step function at T* = 0; 
3) a demonstration that HNC correlations lead to a shell structure (local radial peaks) for 
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strong coupling F > 10 (T* < 0.1); 4) new Monte Carlo simulations to test this theory; 
and 5) a proposal for corrections to the HNC that leads to predictions in good quantitative 
agreement with the Monte Carlo simulations over the range of A^, F tested. 

The qualitative features of the shell structure arising from correlations are: the number 
of shells increases with N, with new shells appearing at special values of A^; a sharpening 
of the shells (more narrow, higher amplitude) with increasing F (decreasing T*); and shell 
populations that grow monotonically with N but which are almost independent of F (T*) (as 



observed in previous MD simulations la, llTI]). The HNC approximation provides agreement 



on all qualitative features of the density profile, e.g. the location of the shells and number 
of particles within the shells is well- described at all coupling values. However, discrepancies 
with the Monte Carlo simulations of the order of 20 — 40% for the widths and amplitudes of 
the shell structure are found. These discrepancies are removed by including a representation 
of the "bridge functions" neglected in the HNC approximation, referred to in the following as 
the adjusted HNC (AHNC). The AHNC theory has the simplicity of the HNC, depending 
only on well-known correlations of the bulk one component plasma, and shows excellent 
agreement with simulations. 

This paper is organized as follows. First the Hamiltonian for the system is defined 
and the density profile is formulated in the Canonical ensemble. The relevant energy and 
length scales are introduced for a dimensionless representation in terms of A^ and F (or 
T*). This is the form most appropriate for Monte Carlo simulation. In Section IIIIl the 
corresponding formulation is given in the Grand Canonical ensemble to allow exploitation of 
density functional theory. A formally exact representation for the density profile is obtained 
in Appendix A in terms of correlations in a bulk one component plasma (OCP), together 
with a corresponding exact equation for the OCP. The HNC approximation is defined in 
Section UTTl as the neglect of the "bridge functions" in these two sets of equations. A further 
neglect of all correlations, the mean field approximation, is explored first over the whole 
temperature range for reference in Section IIVI The effects of correlations within the HNC 
approximation are illustrated in Section |Vl Selected results are compared with those from 
Monte Carlo simulations in Section IVIt where the AHNC is defined and shown to give the 
necessary corrections to HNC required for quantitative agreement. Finally the results are 
summarized and discussed in the last Section. 

In closing this Introduction, it is worth emphasizing that the interesting radial structures 
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studied occur for strong plasma coupling, T > 10. Thus the system presents a good exam- 
ple of a strongly coupled Coulomb system. For very large F the shells become effectively 
systems of particles uniformly distributed on the surfaces of spheres of specific radii. At 
some point the rotational symmetry on these spherical surfaces is broken and these particles 
form spherical Wigner crystals [6]. An adequate description of correlations within density 
functional theory should provide both the fluid phase considered here and the freezing tran- 
sition. The crystal phase is more subtle than in a bulk plasma since it entails packing on a 
non-Euclidean surface, with necessary dischnations depending on the value of considered 
(related to Thomson's problem [fsj]). 

The determination of the non-uniform density for charges confined in a trap is similar 
to the complementary problem of determining the enhanced density in an OCP with an 
impurity ion of the opposite charge. In the latter case the density profile of the bound states 
corresponds to that for the charges in the trap. The formalism here extends to this case in 
a natural way [l^. Both the structure and dynamics for this problem have been studied in 



some detail recently 



191], and will not be considered further here. 



II. MONTE CARLO SIMULATION 

A system of A^ identical charges in a spherical container of radius R is considered. An 
attractive central force is applied at the origin. The Hamiltonian is 

^ /I \ 1 ^ 

H = J2[ ^mvf + Vo in) ) + 2 E ^(^^^•) + 

i=l ^ ^ «5^J=1 

Here Tj and Vj are the position and velocity of charge i. The repulsive interaction potential 
between charges i and j is denoted by V{rij) where Vij = |rj — rj|. The single particle 
potential Vq (^j) represents the central "confinement" potential (referred to as the "trap" 
below) and Vw is the wall potential (zero inside the container and infinite otherwise). The 
focus here is on the average density of charges n{r) for the equilibrium state, with spherical 
symmetry in the fluid phase. In the classical Canonical ensemble this is defined by 

n(ri) = A^^^^- ; — -. (2) 
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Attention will be focused on a harmonic trap and Coulomb interactions among charges. The 
total potential energy in dimensionless form is given by 



r*(rt,..,r^)=/3\/(ri,..,r;v) = r 



,2^3 



N , N 



2g2 2^^^ ^2 ^ |r* -r*| 



V^v^- (3) 



Here, r* = rj/rg, and F is the Coulomb coupling constant F = f3q^/rQ. The usual choice 
for the length scale tq is the ion sphere radius, or mean distance between charges, given 
by 47rrQrl/3 = 1, where n is the spatially averaged density for the confined particles. This 
uniform density can be characterized by the volume of a sphere inside of which the particles 
are localized. This localization can be due to either the walls or the trap alone, depending 
on conditions. Confinement by the trap is determined by the condition that the net force 
on the outermost charge at the surface of this sphere should be zero 

This applies only if Rq < R; otherwise the particles are confined by the wall. In the case of 
interest here, Rq < R and the average density is therefore 



n 



drn{r) = - — -, = -. (5) 



Note that the average density is independent of F and N, although the profile n(r) depends 
strongly on both. 

The dimensionless potential ([3]) is now 

N ^ N ^ ' 

(6) 



and the dimensionless number density is 



N N 



,r- — r ■ 

i=l i¥'j = l I * J 



The parameters characterizing the profile in these units are F and A^. The coupling constant 
F can be interpreted as the inverse temperature in units of the Coulomb energy g^/rg 

The dimensionless density profile n*(r*) can be calculated numerically by direct Monte Carlo 
integration of ([7]) for given F and N. 



- ^ = "BT'i = (8) 



III. THEORY AND APPROXIMATIONS 



For the theoretical analysis of the next sections, it is more convenient to work in the 
grand Canonical ensemble 

^(-0 = ^-'f:W I rfr;..rfr^e-^*(r^'-^), (9) 

iv=o ■ ^ / 
with the normalization constant 

^ = T.^[^) / rfrt..dr^e-^*«-^). (10) 

Here z = e^^, fi is the chemical potential, and A = {h? (3 /2Tim)^^'^ is the thermal wavelength. 
This choice for the ensemble allows application of methods from density functional theory 
The relevant exact equations are summarized in Appendix El In particular, the 



density profile n(r) is related to the excess free energy F^x {P \ n) expressed as a functional 
of that density 

z on [r) 

The excess free energy is a universal functional of the density such that each choice for the 

density corresponds to a different applied external potential. The case of interest here is 

the non-uniform density resulting from the harmonic trap, jSVo (r) —>■ rr*^/2. A different 

case is a uniform density, which results from an applied potential of a uniform neutralizing 



charge density. This latter case is the familiar one component plasma (OCP) [22]. Both are 
described by (11 II) with the same excess free energy functional; only the external potential is 
different. 

In Appendix |X] two exact representations for Fe^ {(3 \ n) are given in terms of a uniform 
reference density. In the first case Vq (r^^) is chosen to be the Coulomb interaction for an ion 
placed at the origin, plus the uniform neutralizing background charge. The density in ( fTTll 
then becomes the pair distribution function 5'ocp(^*) OCP. The appropriate reference 

density is its uniform value far from the ion, for which flTT]) becomes (see Appendix 

^^ghcpin = -rr*-i + J dr*' {ghcpin - 1} cocp (|r* - r*'| ; T) - TBocp {r* \ ghcp) , 

(12) 

together with the Ornstein-Zernicke equation for the direct correlation function cocp 

(ghcpirl - 1) = COCP (r*) + J dr*' {ghcpin - 1} cocp (|r* - r*'| ; T) . (13) 



The function Bqcp {f* I 9ocp) referred to as the OCP "bridge function" [22 1. 

A similar rearrangement of ffTTl) is possible when Vq is chosen to be a harmonic trap using 
a different representation for F^x (/3 | n) , to obtain an exact equation for the trap density 
n(r). In this case the relevant reference density far from the trap center is zero. The analysis 
is outlined in Appendix |X] and the resulting form for Eq. (1111) is 



In- 



-rir*2 + j dv*'n*{r*')cocp (|r* - r*'| ]T)-TB (r* | n*) . (14) 



The second term is expressed in terms of the correlations for the uniform OCP, cqcp i"^*'-, T), 
evaluated at the average trap density n of ([5]). The function B {r* \n*) differs from 
Bqcp ir* I Qocp)^ but by analogy will be referred to as the trap bridge function. Equations 
(fT2|) - f|T4|) provide the formally exact description from which approximations are formulated 
for the trap density. Further details of their context are given in the Appendix. 

Correlations among charges occur in these expressions through the bridge functions and 
cocp- The simplest approximation is the "mean field theory" resulting from the neglect 
of all correlations, B 0,Bocp 0, and cocp{r*) -Tr*~^, for which ([H]) and f|T3]) 
become 

In ^^^^ = -K*^ - J dv*'n*{r*') \t* - r*']'^ , (15) 

r-.te,.M-i)^-.-'-/.."K,.(0-i}k--."r'. (16) 

Equation f|T6|) gives the Debye-Hiickel result for a weakly coupled OCP. The solution to 
( IT5l) for the mean field description of the trap density is discussed below. This mean field 
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14| for the special 



equation was studied for Coulomb and Yukawa traps in references 
case P-^ = T* = 0. 

To account for correlations the HNC approximation is used, which results from neglecting 
the bridge functions B 0,Bocp but retaining the effects of cqcp (which in this 
approximation is denoted chnc) 



In = -P-r*2 + j dr*'n*{r*')cHNC (|r* - r*'| ; P) , (17) 

ln^^^^(r*) = -Pr*-i + J dr*' - 1} Chnc (|r* - r*'| ; P) , (18) 

idHNcirl - 1) = CHNC in + I dr*' {g*HNcin - 1} CHNC (|r* - r*'| ; P) . (19) 
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Equations (fT8|) and (fTOjl are a closed set of equations to determine gnNC chnc for the 



OCP, a well-studied problem [2^. Equation f|T7|) is an extension of the HNC to an approx- 
imation for the localized trap density. It is shown below that the inclusion of correlations 
leads to qualitatively different results for the density compared to results from the mean 
field theory. 

To go beyond the HNC approximation requires some estimate of the bridge functions. 
This is a difficult problem in general. For the OCP, effects of Bqcp are important only at 
very strong coupling. However, for the trap density the bridge functions have a quantitative 
importance even at moderate coupling. This is discussed in Section IVTl where a phenomeno- 
logical estimate is explored and shown to remove the observed discrepancies between HNC 
and Monte Carlo simulation results. Further work in this direction is planned for future 
studies. 

In the following it is convenient to write the above equations using a representation for 

the density in terms of a dimensionless effective potential U* (r*), defined by 

_ g-rt/*(r-*) 
n*ir*) = N . (20) 

Here the chemical potential /x has been eliminated in favor of the average number of particles 
N = f drnir). Substitution of ( l20l) into (fT4l) gives the corresponding exact equation for 

U* (r*) 

1 ^ dr*'e~^^'^''*''^^)ci\v* -r*'\ Fl 
U*{r\V,N)=-r*^ + N^-^^ ' +B{r\n), (21) 

^ ' ' ^ 2 j d^*,^-^U*{r*'V,N) \ \ >^ \ } 

with the notation c(r*, F) = — Fc(r*, F). The solution is parameterized by only two dimen- 
sionless constants, F, A^. 



IV. MEAN FIELD THEORY 



To develop a quahtative understanding of the dependence of the profile on the parameters 
F and it is useful to explore this simple mean field approximation in the dimensionless 
variables of the last section. Taking the mean field limit of (12T|) . and performing the angular 
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integration in the second term gives the mean field equation 

1 TV 

U* {r*,T,N) = -r*2 + 



2 /;^*dr*r*2e-r^*r'^'^ 

r* R* 

I J dr*' r-^e-r^l'-'.r.^) + J dr*' r*'e-^^*(^*'^-^) J . (22) 



X 







For high temperatures T* = T ^ >> 1 (l22l) gives 

f/.(..,r.]v)^i(i-^).- + §, (23) 

where an outer hard wall at R* has been restored. Consequently the trap density becomes 

/ — V N pi /'i ^ V*2 

n*(r*,r,N)^ J — ^^e'^^y^'T^r . 24 



It is understood that r* < R* and the density is zero otherwise. For large R* fl24j) is just 
the expected Boltzmann factor for uncorrelated particles in an external harmonic potential. 
However, as the number of particles in the trap increases their Coulomb repulsion competes 
with the confining trap potential, enhancing the density toward the outer wall. For < R*^ 
the density increases toward the center of the trap, while for > R*^ it increases in the 
direction of the outer wall. The density is uniform at the special value R* = N = Rq, the 
point where the Coulomb repulsive force and attractive confinement force exactly balance. 

At low temperatures T* = T^^ « 1 and e"'"'^*^'''-* — in (l22l) unless U* (r*) is constant, 
in which case these factors cancel in the numerator and denominator. Therefore, a solution 
is sought of the form 

{II* r* <^ r* 
0, max ^ ^25) 

U* (r*) > 0, r*> C 

for some constant and r*,.,„ < R*. Then, for F >> 1 fl22D gives for r* < r* „^ 



[/•(r-.r,¥)^i(l--^)r- + i^ (26) 



max 



The first term must vanish for consistency with ( !25l) . which determines r^g^-^ 

C = (27) 
Next, for F » 1 and R*o < r* < R* ^ gives 

U*{r*,T,N)^^-r*' + ^. (28) 
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FIG. 1: (color online) Mean field results for = 100 for various values of T. The Coulomb limit 
(r — > oo) is also shown. 



Consequently, the trap density in this limit is 



n (r 



3N 



_r( ir*2 + JV 

e 



r* > R*Q 



(29) 



The density is uniform up to -Rq which depends only on A^, and is exponentially small for 
larger r*; the density effectively has a step profile. 

More generally, (!22l) can be solved numerically for the full range of T and A^. The 
dependence on F is illustrated in Figure 1 for = 100 and F = 0.1, 1, 10, and 20. Also 
shown is the limiting step function for F ^ oo. The edge of the profile is the radius at which 

1/3 

the confinement force equals the total Coulomb force, Rq = N . Although instructive as 
a reference, the low temperature behavior within the mean field approximation is relevant 
only for the average density, since charge correlation effects become dominant for F > 10, 
as shown in the next section. 
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I-H 




FIG. 2: (color online) Dependence of the scaled direct correlation function on T. Also shown is the 
small r Coulomb limit (mean- field approximation ). 

V. HNC APPROXIMATION 

The density in the HNC approximation is given by (|T7|) - (|T9l) . The effective potential is 
determined from the equation 



U* {r*,T,N) 



Jdr 



*f-ru*(r*',r,N' 



CHNc{\r* 



r*'|,r) 



.ru*{r*'r,N) 



(30) 



J dr*'e 

The only difference from the mean field form is the replacement of the Coulomb interaction 
by the direct correlation function chnc- This is first calculated from the corresponding HNC 
approximation for the OCP, ( |T8i) and ( fT9i) . and the results are illustrated in Figure 2 for 
F = 0.1, 1, 10, and 100. It is seen that the Coulomb limit is recovered for all r* in the limit 
F ^ 0, and for any F at sufficiently large r*. 

With Chnc known, fl30|) can be solved to determine the trap density. Figure 3 shows the 
results for the same conditions as in Figure 1, illustrating the effects of correlations relative 
to the mean field results. The monotonic decrease with r* is now replaced by the appearance 
of structure ("shells") for F = 10 and 20. The latter two curves also suggest that the location 
of the shells is insensitive to the value of F, as is confirmed below. Figure 4, for F = 40, 
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FIG. 3: (color online) Radial density profile with correlations included, for the same conditions as 
Figure 1. 
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FIG. 5: (color online) Shell population for two-shell conditions. The higher trend curve is for the 
the outer shell. Data is shown for HNC and MC calculations. 



shows that the number of "shells" increases with A^. Figure 5 illustrates the formation and 
filling of the first three shells, giving the number of particles in each shell (defined as the 
number between two spheres at successive radial minima) as a function of for F = 10, 20, 
and 40. The initial dependence is linear as only a single outer shell is populated. However, 
after formation of the second shell the dependence is more complex as additional particles 
can go into either shell. Remarkably, this process of filling is almost independent of F (in 



agreement with experiment and simulations 
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VI. COMPARISON TO MONTE CARLO SIMULATION 

To assess the quality of predictions based on the HNC approximation some benchmark 
Monte Carlo simulations have been performed for comparison. For moderate coupling, 
F = 1, correlations are important (see Figure 2, for example) but shell structure does 
not appear. Figure 6 shows the comparison with Monte Carlo results for F = 1, and 

= 10, 100, and 1000. Not shown are results for < 10 where the agreement is poor even 
at weak coupling. The first conclusion is that the HNC approximation is accurate at weak 
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FIG. 6: (color online) Comparison with Monte Carlo results for F = 1. 

to moderate coupling, F < 1, and > 10. 

The formation of shell structure at larger values of F is tested in Figure 7 for F = 20, 
= 25 and 100, and in Figure 8 for F = 40, = 300. It is seen that the formation and 

location of the peaks is very well described, but systematically their widths are too large 

and amplitudes too small by as much as 20 — 40%. Nevertheless, the particle number in the 

shells is quite accurate, as shown in Figure [51 



VII. ADJUSTED HNC 



The HNC approximation provides a surprisingly good representation of correlation effects 
on the density profile of strongly coupled charges in a trap: formation, population, location, 
and number of shells. Its primary limitation is the significant discrepancy in the amplitudes 
and widths of the shells. This is due to neglect of important contributions from the bridge 
functions B {r \ n) and Bocp{r \ in ( JT2i) and ffl31l at strong coupling. This problem 
was addressed some time ago for the OCP by Ng [231], who observed that the peaks of the 
HNC pair distribution function gocpi^*) appear at the correct position but with incorrect 
amplitudes and widths. This is clearly analogous to the current errors of the trap density. 
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FIG. 7: (color online) Comparison with Monte Carlo results for T = 20. 

Ng corrected the HNC by an estimate of the OCP bridge function in the form 

Bocp (r I n) ^ A (F) K. (r*) , K) (r*) = ^, (31) 

where A (F) is a chosen function to interpolate between A (0) = and some constant value 
A {oo). The advantage of the chosen form is that when inserted in f[T^ the form of the HNC 
approximation is recovered, but with a renormalized coupling constant F' = (1 + A (F)) F 

lng*ocpin = -FV*-i + J dr*' {ghcpir*') - 1) Cocp (|r* - r*'| ; F') . (32) 

Note that Cqcp ^ ; F') also becomes a function of the renormalized coupling constant 
since it inherits that dependence only from 5'ocp('"*) Ornstein-Zernicke equation (fT3l) . 

Ng found that the Monte Carlo data for gocpi^*) could be fit within two percent with the 
choice A (oo) = 0.6. Alternatively, (131 p can be viewed as a representation of the OCP bridge 
function determined from f[T^ using Qocpi^*) from Monte Carlo as input. 

It is plausible to use this same scheme for the trap bridge function B {r \ n) as well 

B{r\n)^X (F) Vo (r*) , (r*) = ^r*\ (33) 

with the same choice for A (F). Accordingly, f|T^ gives the HNC form for the trap density 

In ^^^^ = -F'-r*2 + f dr*'n*{r*')cocp (|r* - r*'| ; F') . (34) 
z 2, J 
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FIG. 8: (color online) Comparison with Monte Carlo results for N = 300 and T = 40. 

This is the same as the HNC approximation of Section |Vl except with the bridge function 
corrected value T'. Based on the resuhs of Ng, the same value of A (oo) = 0.6 for the 
strong coupling conditions is used here. Equations fl32l) and flM|) will be referred to as the 
adjusted HNC (AHNC) approximation. Figures 9 and 10 show a revised comparison for the 
conditions of figures 7 and 8 but now with the AHNC calculated with F' = 1.6F. Most of the 
discrepancies are seen to be removed by this AHNC approximation. It describes accurately 
the complexities of the trap density over a range of F, including strong coupling, with 
the structural simplicity of the HNC approximation and only OCP correlations as input. 
Its domain of applicability and limitations will be explored in more detail elsewhere. It is 
somewhat remarkable that the simple forms fl3Tl) and fl33l) for the bridge functions work 
so well. This success should motivate future work to address the theoretical basis for such 
practical approximations to the bridge functions. 

VIII. DISCUSSION 

The objective here has been the development and exploration of a theoretical framework 
for analysis of the density profile of charges in a harmonic trap. Previous extensive work 
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FIG. 9: (color online) Comparison of adjusted HNC with Monte Carlo results for F = 20. 




01 2345678 

r* 



FIG. 10: (color online) Comparison of adjusted HNC with Monte Carlo results for N = 300 and 
F = 40. 
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on this problem has focused on experimental and simulation studies of the low temperature 
crystal structure (and its melting) for Coulomb, Yukawa, and Lennard- Jones interactions. 
Here the complementary theoretical study starting from the high temperature fluid phase 
has been described, showing the evolution of shell structure at lower temperatures due 
to correlations that are precursors of the observed structure in the crystal phase. A first 
order representation of correlations in the HNC approximation is shown to capture the 
radial structural features of the density profile. Particularly remarkable is the formation 
of "blurred" shells for moderate coupling at the magic numbers for tiling the sharp crystal 
shells. This anticipation in the fluid phase of the low temperature crystal structure will be 
discussed elsewhere. 

The analysis above provides justification for a simple and practical description within 
the fluid phase for charged particle confinement, the AHNC approximation. Although the 
discussion here has been limited to the case of Coulomb interactions there is no complication 
involved in extending it to other more practical forms such as Yukawa particles. The formu- 
lation of this approximation within the context of density functional theory also provides a 
complete thermodynamics for this system. In this way it is hoped that the AHNC could be 
improved to include a description of freezing and crystal structure. 
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APPENDIX A: REPRESENTATIONS FROM DENSITY FUNCTIONAL THE- 
ORY 

Consider a given system in an external potential. For simplicity only the case of spher- 
ically symmetric central potentials are considered. Two classes of external potential are 
identified. In the first class, the potential goes to zero or a constant at large r, resulting in 
an average density that becomes uniform at large distances. This is the usual case of local- 
ized perturbations in a bulk fluid. The second class is potentials that become unbounded 
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for large r, so that the resulting average density is localized in some region about the ori- 
gin. This is the case of interest here for charges in a trap. In the following, different exact 
representations of the inhomogeneous density appropriate for each of these cases are given. 



The approach is based on the formal structure of density functional theory 
review briefly the relevant equations, the grand potential is defined first by 



2l|. To 



pn{p\ fx-Vo) = -\nEj2 Tre-^^^-^'^\ (Al) 

where S is the grand partition function of ffTOl) . The notation makes explicit the functional 
dependence of Vt on the external potential Vq. The density is obtained from the functional 
derivative of VL 

The roles of n{r) and Vo(r) can be exchanged by the Legendre transformation 

F{P\n)=n{(3\^i-V^) + J dr{fi-Vo{r))n{r), (A3) 

The new functional F {(3 \ n) is the free energy expressed as a universal functional of the 
density. Equations flA2p and flA4p make explicit the fact that the density is a functional of 
the applied potential and vice versa. The functional F {(3 \ n) is universal in the sense that 
each choice for Vq corresponds to evaluation of this functional at a corresponding specific 
density field. Equation ( ]A4I) will be applied here as an exact equation for the density in 
terms of the chosen Vq. 

The free energy can be divided into an "ideal gas" form Fo{l3 \ n), (its form in the absence 
of interactions among the particles), plus the "excess free energy" due to interactions, F^xiP \ 
n) 

F{f]\n) = Fo{f3\n) + F,,{(3\n). (A5) 
The first term can be calculated directly so that (1A4I) takes the more practical form of (ITTll 



z dn[r) 

A further decomposition of F^.^ is useful for the introduction of approximations. The 
central idea is to exploit knowledge of correlations in uniform systems to characterize the 
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non- uniform density profile due to Vq. First, note tlie identity for any functional X{g) of 
the function g (r) 

X{g)-X{gr)= [' dXd^X [Xg + [1 - X) gr] 
Jo 

= dX j dr^-^^ |Ag+(i-A)g, (r) - (r)) (AT) 

where g^ (r) is an arbitrary reference function. Application of this identity twice to the 
excess free energy F^x gives 



\nr] n r — rir 



F,x{f3 I n) = F,^{(3 \nr)- j drp-'c^'^ (r 

- [ dX [ dX' [ drdr'p'^c^^^ [r, r'\X'n + (1 - A') n,] (n (r) - n,) (n (r') - n,) . 
Jo Jo J 

(A8) 

Here rir is an arbitrary reference density, and the direct correlation functions have been 
introduced by the definitions 

& ^ ri, r2..r„ n) = -(3 — -. (A9) 

dn {ri)dn (ra) ..dn (r„) 

For external potentials Vq that vanish for large r there is a natural reference density n 
corresponding to the uniform system at distances far from the perturbation. Then flA2p 
becomes 

F,,{P\n)=F,,{P\n)~ [ rfA (1 - A) 

X J dvdr'p-^c^^^ [r, r'lA'n + (1 - A') n] (n (r) - n) {n (r') - n) . (AlO) 



Here the second term of ( 1A8I) vanishes by requiring that nV = N, the average total particle 
number. Again applying (1A7I) to express c^^-* (r, r'|A'?7. + {1 — X')n) in terms of the uniform 
density gives 

F,x{f3 I n) = F,x{f3 drdr'p-'c^^^ (|r - r'| \n) {n (r) - n) {n (r') - n) 

-I dX{l-X) I dX' [ drdr'dr" (n (r) - n) {n (r') - n) {n (r") - n) 
Jo Jo J 

X c(^) [r, r', r" \ X'n + (1 - X')n] (All) 
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Next, consider external potentials that become unbounded for large r so that the density 
vanishes at large distances. In this case the appropriate choice for the reference density in 
[8]) is = 0, giving 

Fe^{p\n) = dX{l-X) J drdv'n (r) n (r') /J-^c^^) (r, r' | An) , (A12) 



The first two terms on the right side of (lASP vanish at zero density, and the first A integral has 
been performed. The integration domains are now controlled by the density, which restricts 
them to the bounded domain. In this domain, it is reasonable to identify an appropriate 
uniform reference density n and to express c*^^^ (r, r' | An) in terms of this uniform density 
by a further application of (1A7I) 



c*-^-* (r, r' I An) = c*-^^ (r, r' | n) + / dxd^c^'^^ [r, r' | xXn + (1 — x)n] 

Jo 



c*-^-* (r, r' I n) 



+ dx j dr" (An (r") - n) c^^^ [r, r', r" | a;An + (1 - x)n] . (A13) 
Then the excess free energy becomes 

FexiP \^) = -\j d'^d'^'n (r) n (r) /J^^c^^^ (|r - r'| | n) 

- f dX{l-X) f dx f dr"p~^c p^r, r', r" | xAn + (1 - x)n\ ((An (r") - n)) , 
Jo Jo J 

(A14) 

Equations ( lAlip and ( ]A14I) are both formally exact and equivalent, but each is formulated 
so that the leading terms are reasonable approximations for different classes of external 
potentials. Use of flATTj) in (jA6l) leads to the form ([T2D of the text; use of flATil) in (lA6|l 
leads to the form f|T^ . The Ornstein-Zernicke equation ([T3|) is an identity that follows from 
the definition of c^^-* in (]A9p . The last terms of (lAlip and flA14p give rise to what is referred 
to as bridge functions in the text, and provide the formal definitions for each. 

The above is quite general. For the case of interest here, a system of charges, the uniform 
density case corresponds to the one component plasma. 
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